Continuum Electromechanical Modeling of Protein-Membrane Interaction 
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A continuum electromechanical model is proposed to describe the membrane curvature induced by 
electrostatic interactions in a solvated protein-membrane system. The model couples the macroscopic 
strain energy of membrane and the electrostatic solvation energy of the system, and equilibrium membrane 
deformation is obtained by minimizing the electro-elastic energy functional with respect to the dielectric 
interface. The model is illustrated with the systems with increasing geometry complexity and captures the 
sensitivity of membrane curvature to the permanent and mobile charge distributions. 
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Boundaries of eukaryotic cells and most organelles 
within the cells are defined by lipid bilayer membranes. 
Shape and topological transformations of membrane are 
crucial steps in numerous transport and signaling processes 
of cells, including cell migration, membrane trafficking, 
and ion conductance ll|-|3D. There are various types of 
proteins that are displaced to the vicinity of membranes 
or are embedded in them. Many of the delicate changes 
of configuration and topology of biological membranes re- 
sult from the interactions between lipids as well as be- 
tween lipids and proteins. Among the most interesting ex- 
amples of protein-induced membrane deformation are the 
Bin/amphiphysin/Rvs (BAR) domain-induced membrane 
curvature [4] and the Endosomal Sorting Complex Re- 
quired for Transport III (ESCRT-III) induced membrane 
budding or protrusion These interactions happen over a 
wide range of time and length scales, and depend critically 
on the lipid composition, charge distribution, hydrophobic - 
ity, and elastic moduli. Understanding of quantitative rela- 
tionships between the forces and the membrane geometry 
have long been subjected to intensive studies using either 
discrete methods such as molecular dynamics (MD) simu- 
lations and the coarsed grained MD ]6|1, or continuum elas- 
ticity based on minimization of Helfrich-Canham-Evans 
bending free energy and its simplifications [9], some- 
times under the constraints such as surface area or enclosed 
volume ifToll . While the dynamics in the lipid-protein inter- 
actions can be revealed in atomistic details by MD simu- 
lations, protein-induced macroscopic bending, expansion, 
contraction of lipid bilayers can be more efficiently simu- 
lated using continuum mechanics at substantially reduced 
computational cost. The continuum models also give us 
access to interacting at the interfaces between molecular 
biology, mathematical modeling, and numerical comput- 
ing. Most of the current studies of membrane mechanics 



model the lipid bilayer as an elastic sheet with vanishing 
thickness, and are focused on the the membrane's curva- 
ture energy, stability, or the deformation modes for given 
external mechanical constraints or loads. In this Letter, we 
propose a continuum electromechanical energy and asso- 
ciated system of nonlinear partial differential equations as 
a framework for modeling the self-consistent macromolec- 
ular deformation induced by the long range electrostatic 
interactions in solvated macromolecular systems. 

The essential idea underlying our model is the coupling 
of the electrostatic potential energy of the solvated macro- 
molecular system and the nonlinear elastic energy of flexi- 
ble biomolecules modeled as elastic continuum. We stress 
that, as opposed to other continuum models, the bilayer is 
described as an elastic sheet with a finite thickness and a 
constant dielectric permittivity. We neglect the atomistic 
details of the lipids, but the partial charges in the head- 
groups of the lipid bilayers, modeled as a distribution of 
charge on the surfaces, will be retained and incorporated 
into the modeling of the electrostatic interactions with the 
proteins, c.f. Fig. [T] In contrast, the proteins are modeled 
as rigid bodies with atomistic detail. Their dielectric in- 
terfaces with the solvent are described by properly defined 
molecular surface. 

Let (j) be the electrostatic potential of the entire solvated 
system and u be the displacement vector of its flexible sub- 
domain with a mass density p m , the electro-elastic energy 
density of the system is defined to be 
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where the first and the second terms are the kinetic and po- 
tential energy of the elastic domain in the system, respec- 
tively, and e = |(Vii + Vu T ) is the linear strain tensor. 
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FIG. 1. Illustration of the proposed electromechanical model 
for protein-membrane interactions. Left: Molecular surface of 
a protein-membrane complex. Atomistic details are retained in 
the protein but are neglected in membrane. The left and right 
surfaces of the membrane are negatively charged, and its circu- 
lar face is charge free and constrained. Depending on the mobile 
ion concentration and the partial charge distribution in the pro- 
tein, the membrane might be attracted to bend towards the protein 
(Middle) or repulsed away from the protein (Right). 

A, [i are the two Lame constants. The electrostatic poten- 
tial energy (density) is given by g, for which we use the 
following ansatz: 
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where e is the spatial-dependent dielectric permittivity, c$ 
is the concentration of the zth species of mobile ions in the 
solvent with charge The permanent charge distribution 
p q consists of the M singular point charges located at x t in 
the proteins and the constant charge distribution p s on the 
surfaces S m of the membrane, i.e., p q = J2i=i Qi^( x i) + 
p s [x(S m )], where 5 and \ are me 3-D and 2-D Dirac delta 
functions, respectively. It is known ifllll that this potential 
(j> can be solved from the Poisson-Boltzmann (PB) equa- 
tion, -V • (eV0) - Y!Li c i q l e- q ^ /kBT = p q , and the 
electrostatic force density / = — Vg is given by 
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where E = — V0, T is the molecular surface, and the sub- 
scripts s, m denote the values in the solvent and in the 
molecules, respectively. Here we choose e s = 80 and 
e m = 1. The surface charge distribution will be treated 
as an interface condition on the membrane surfaces, i.e., 
— (e s V</> s — e m V0 m ) • n = p s , where n is the surface 
outer normal. We finally consider the energy functional 



S= / Fdxdt (4) 
Jti Jn 

and its variation with respect to the spatial displacement 
5S = J^ 2 J n STdxdt. The minimization of the energy 
functional gives rise to an equation 8 J- = 0, i.e., 
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where I is the identity matrix. This shows that the min- 
imization of the electro-elastic energy functional is 
equivalent to the coupled solution of the PB equation and 
the elastic equation (01. The coupling is realized through 
supplying the electrostatic force computed from the PB 
equation to the elastic equation as a dynamical load and 
supplying the varying membrane configuration computed 
from elastic equation to the PB equation as a dielectric in- 
terface. The time derivative in the elastic equation can be 
removed if one is only interested in the equilibrium con- 
figuration of the deformable macromolecules. This regen- 
erates the direct coupling of the nonlinear elasticity equa- 
tion and the PB equation 11211 . and gives us a handle of 
high-fidelity numerical methods for solving these partial 
differential equations. Furthermore, the equivalent relation 
between the electrostatic force density / and the Maxwell 
stress tensor (MST) [13] allows us to use the MST to com- 
pute the electrostatic force exerted only on the surfaces of 
the membrane: 
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f = M-n- k B Tj2ci(e~ 9 ' HkBT ~ 1), 

M = -|e s |V0 s | 2 + e s V4> s <g> V0 S . 

It is worth noting that the charged membrane will subject 
to a self-force / in the absence of the proteins. This self- 
force must be subtracted and thus the net traction applied 
on the membrane is / — f . 

To ensure the accuracy of the numerical solution of elec- 
trostatic potential induced by singular point charges, sin- 
gular surface charge distribution, and mobile ions, we in- 
troduce a stable decomposition of the PB equation 111411 . 
with which one can analytically solve the potential compo- 
nent induced by the singular charges. A molecular surface 
conforming finite element method is then used to solve the 
regular component of the potential with which the interface 
conditions due to the surface charge distribution can be en- 
forced exactly. In order to avoid regenerating interface- 
conforming finite element meshes due to the membrane 
deformation in the iterations, a Piola transformation Ill5ll 
is adopted so that the elastic equation can be solved on the 
same reference configuration with varying surface traction. 

As a proof of the model and the numerical method, we 
present a test system in which a cylindrical membrane is 
buckled by a model protein whose structure is adjustable. 
The membrane has a radius of 150A and a thickness of 
30A. The radii of atoms in the model protein are uniformly 
lOA. The default surface charge density — 0.0044e/A 2 is 
computed by assuming that the membrane has the same 
composition of 30% dioleoylphosphatidylserine (DOPS) 
and 70% dioleoylphosphatidylcholine (DOPC) as studied 
in lloll lflr^l . We fix the Young's modulus and the Poisson 
ratio of the membrane to be lpN/A 2 and 0.3, respectively 
[fTvh . The number, positions and charges of the atoms in 
the protein, the membrane surface charge density, and the 
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mobile ion concentrations will be varied to validate the pro- 
posed model for simulating the membrane deformation un- 
der a wide range of electrostatic forces. Fig. [2] plots the 
membrane displacement as a function of stated parameters. 
The plots show that rich dynamics can take place in this 
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FIG. 2. Change of membrane displacement with different pa- 
rameters for the test systems. Unless specified particularly the 
charge of atom is 5e, the mobile ion concentration is 150mM, and 
the surface charge density assumes the default value. In the first 
system (A,B,C) the protein consists of a single atom centered at 
(0, 0, 50). In the second problem (D) the overall electrically neu- 
tral protein consists of 8 model atoms centered at (±6, ±6, 50) 
and (±6, ±6, 56). Every atom on z — 50 and z — 56 carries a 
charge of 2 and —2, respectively. The center of the membrane is 
initially placed at (0, 0, 0) in both problems with its axle aligned 
with the z-axis. 

electromechanical system and different types of electro- 
static force compete. A large negative surface charge den- 
sity results in a strong attractive interaction with the pos- 
itively charged model protein, and gives rise to a positive 
displacement. In the absence of membrane surface charge, 
a small but negative displacement (Fig. |2j\) suggests that 
the total electrostatic force induced by the discontinuous 
dielectric permittivity and the ionic osmotic pressure is re- 
pulsive on the low dielectric material in agreement with the 
analysis in jl2lll8ll . The attractive interaction gets stronger 
with the increase of the net positive charge of the protein, 
resulting in a monotonical increase of the membrane dis- 
placement (Fig. |2j3). The mobile ion concentration has a 
more complicated consequence in regulating the protein- 
membrane electrostatic interaction (Fig. WP)- The present 
of mobile ion at low concentration reduces the electrostatic 
potential and in turn the electrostatic force on the mem- 
brane surfaces, resulting a decease of the displacement. 
Nevertheless, the further increase of the mobile ion concen- 
tration from 0.8mM to 2.2mM leads to a slight increase 
of the membrane displacement. This suggests that over this 
small range of mobile ion concentration the decrease of 
the attractive force is dominated by the decrease of repul- 
sive force components, resulting in an increase of attractive 
net force. This dominance reverses with the further im- 



provement of the ion concentration, as seen in the decrease 
of the displacement for all ion concentrations larger than 
2.2mM. In the second test system, the lower positively 
charged domain of the protein has an attractive interaction 
stronger than the repulsive interaction between the nega- 
tively charged upper domain and the membrane. One the 
other hand, the repulsive interaction due to the discontinu- 
ity of the dielectric permittivity is also prominent under the 
weak screening of the mobile ions, and thus a negative dis- 
placement as large as — 8 A is seen in Fig.|2jD. Noticing that 
the dielectric pressure (2nd term in Eq.©) is proportional 
to the square of the magnitude of the electric field while 
the qE force (1st term in Eq.©) is proportional to E, this 
pressure reduces more quickly as the electric field weakens 
with the increase of ionic strength. At a sufficiently high 
ion concentration close to the physiological condition, a 
small positive displacement is found. 

We now investigate the membrane curvature induced by 
the BAR-domain protein (PDB ID: 1I4D) dimer and com- 
pare our results to the molecular dynamics simulations by 
Blood et. al. [6]. The membrane at free state is a rectangu- 
lar cuboid whose size is (x x y x z) = [—100, 100] x 
[—50,50] x [— 15,15]A 3 . The concave surface of the 
dimer has a radius of curvature 80A. The coordinates of 
the dimer structure are rotated and translated such that the 
three atoms (GLN225:0,LEU149:HG,SER164:HG) are on 
the x — z plane and the z -coordinates of the last two atoms 
are 0. The z-distance between these two atoms and the top 
surface of the membrane, which is initially placed below 
the dimer, is denoted by D. The top and bottom membrane 
faces are charged at the default value, and the other four 
faces are charge free. The two end faces in ^-direction are 
fixed and thus homogeneous Dirichlet boundary conditions 
will apply. The periodic boundary condition is enforced in 
y-direction to reduce the numerical artifacts due to finite 
truncation. The results for the membrane displacement and 
curvature, measured along the intersection between the top 
surface and the plane y = 0, are shown in Fig.[5J The dimer 
has a positively charged concave face and a net charge of 
— lOe. This negative net charge gives rise to a strong repul- 
sive qE force to the negatively charged membrane surfaces 
when the ion screening is weak, which along with the re- 
pulsive dielectric pressure leads to a positive bending cur- 
vature of 2.2 x 10~ 3 /A. A transitional ion concentration is 
found at 8mM, below which the bending curvature is pos- 
itive while above which the curvature is negative. Similar 
to the second model problem we tested above, this transi- 
tional ion concentration corresponds to the state in which 
the screened attractive interaction between the positively 
charged side of the dimer and the negatively charged mem- 
brane surfaces is balanced by the repulsive interaction be- 
tween the negative charged side of the dimer and the mem- 
brane surfaces plus the dielectric pressure as well as the 
ionic pressure. The maximum curvature, 3.7 x 10~ 3 /A, 
occurs at an ion concentration of 200mM, while the cur- 
vature at the physiological concentration, lOOmM in this 
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FIG. 3. First row: The structure of protein 1I4D (Left); Sur- 
face triangulation of the BAR-domain protein-membrane com- 
plex (Middle) and the electrostatic potential on the protein and 
membrane surfaces (Right). Second row: Membrane curvature at 
various ion concentrations (Left) and membrane displacement at 
various initial separations D (Right). 

study, is 3.2 x 10 _3 /A. These demonstrate that the present 
model is able to capture the critical role of the mobile ions 
in mediating the protein-membrane interactions. Neverthe- 
less, all these bending curvatures are smaller than the in- 
trinsic curvature of the dimer (1.25 x 10~ 2 /A). We specu- 
late that a larger curvature could be found at suitable initial 
separation between the membrane and the dimer. However, 
the simulations with various separations (c.f. Fig. [5]) illus- 
trate that the bending curvature becomes slightly smaller 
as the separation gets smaller. These observations sug- 
gest that the attractive electrostatic force is not sufficiently 
strong to generate a bending curvature as large as the in- 
trinsic curvature of the dimer when the membrane is con- 
strained at two a>ends. On the other hand, a charged mem- 
brane free of these constraints will move towards the dimer 
under the attraction and eventually contact with the two 
ends of the dimer; these two ends will then serve as the 
constraints for the bending of the membrane. This pro- 
cedure has been observed in molecular dynamical simula- 
tions Jfl @], where the two ends of the dimer will insert 
into the membrane after contacting with the membrane and 
further facilitate the bending of the membrane. 

The results clearly demonstrate that our continuum elec- 
tromechanical model efficiently captures the balance be- 
tween the electrostatic energy and the elastic strain energy 
in the protein-membrane interaction. The model quanti- 
tatively characterizes the dependence of the interaction on 
the solvation, partial charge distribution in protein, mem- 
brane composition and its surface charge distribution, and 
the membrane dielectric constant, some of which has been 
observed in previous studies lalSD, but never described in 
a continuum framework. The present formalism is mainly 
limited by the fixed constraints on some un-charged mem- 
brane surfaces, which could be resolved by introducing the 



protein surface as a limit of membrane deformation ifl^l . 
Future efforts to improve these approximations are critical 
to accurately describe the protein-membrane interactions 
for more complicated geometries. The formulation and ex- 
amples presented here are expected to spur the interests of 
other investigators in molecular biology, mechanics, and 
mathematics. 
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